home *** CD-ROM | disk | FTP | other *** search
- {*********************************************************************
- * 円周率の計算 *
- * HAPPyのサンプルプログラム (TOWNSフリコレ6版) *
- * Auther H.Asano *
- *********************************************************************}
-
- { 定数AREA-7桁分の円周率を求めます。
- しかし、不思議なことに、小数点以下11桁目の答えが違うのです。
- 正解は8なのに、このプログラムでは9となります。
- 1000桁まで計算すると、761,762桁目も違うことが判明しました。
- どなたか、このプログラムの誤りを指摘して下さい }
-
- { マーチンの公式 π/4 := 4*arctan(1/5) - arctan(1/239) を 使う。
- π := 16*arctan(1/5) - 4*arctan(1/239)
- arctan(x) := x/1 - x*x*x/3 + x*x*x*x*x/5 - ・・・
- x/1の分子分母にxをかけて整理すると
- 16*arctan(1/5) := 16*5/25/1 - (term)/25/3 + (term)/25/5 - ・・・・
- = 80 /25/1 - ・・・
- -4*arctan(1/239) := -4*239/57121/1 + (term)/57121/3 - (term)/57121/5 + ・・・
- = -956 /57121/1 + ・・・
- ここで、termは、前項までのx指数計算の結果を格納している。
- たとえば第1項の80/25が次の項のtermである。第2項は 80/25/25 / 3 となる。
- }
-
- program Pai(output) ;
-
- const AREA = 187 ; { 180桁まで求める }
- {***** このプログラムをHAPPyのスピードアップ測定に使ってきました。
- 使用マシンは、FM TOWNS 40Hです。
- 以下は、スピードアップの履歴です *****}
- { 1992.4.13 1分55秒 }
- { 1992.4.18 1分49秒 }
- { 1992.4.18 通訳系の最適化により 1分44秒 }
- { 1992.4.22 通訳系の書き直しにより 1分14秒 }
- { 1992.4.30 直訳系と通訳系の分離 1分10秒 }
- { 1992.5.01 通訳系の記憶装置型変更 1分00秒 }
- { 1992.5.02 通訳系の改良 0分58秒 }
- { 1992.5.02 program の改良 0分56秒 }
- { 1992.11.03 通訳系の改良 0分54秒 }
- { 1992.11.26 通訳系の改良 0分48秒 }
- { 1993.02.17 programの改良 0分38秒 }
-
- type range = 0..AREA ;
- wkplace = array[range] of integer ; { 添字2と3の間が小数点 }
-
- var pai : wkplace ; { 円周率の格納エリア }
- work,term : wkplace ; { 計算上のワーク }
- i : range ;
- p : range ; { 収束位置 }
-
- {********** 割り算処理 **********}
- { wk := term div jyosu を 計算する。
- pは、termi配列の0でない最初の位置を示す。これは収束判定に使う }
- procedure Divide(var wk : wkplace ; jyosu : integer) ;
- var i : range ;
- hi,w : integer ;
- begin
- { 有効数字が出てくる位置を求める }
- while (term[p] = 0) and (p < AREA) do p := succ(p) ;
-
- hi := 0 ;
- for i := p to AREA do
- begin
- hi := hi*10 + term[i] ;
- w := hi div jyosu ;
- if w >= 1 then
- begin wk[i] := w ;
- hi := hi mod jyosu
- end
- else wk[i] := 0
- end
- end ;
-
- {********** 加減算 **********}
- { pai := pai + work ・・・・ AddSubFlag = true の時
- pai := pai - work ・・・・ AddSubFlag = false の時 }
- procedure AddSub(AddSubFlag : Boolean) ;
- var carry : Boolean ; { 次の桁に影響する時 true }
- wk : integer ;
- i : range ;
- begin
- carry := false ;
-
- if AddSubFlag then
- begin { 加算 }
- for i:=AREA downto p do
- begin
- wk := pai[i] + work[i] + ord(carry) ;
- carry := wk >= 10 ;
- if carry then pai[i] := wk - 10
- else pai[i] := wk
- end
- end
-
- else
- begin { 減算 }
- for i:=AREA downto p do
- begin
- wk := pai[i] - (work[i] + ord(carry)) ;
- carry := wk < 0 ;
- if carry then pai[i] := wk + 10
- else pai[i] := wk
- end
- end
- end ;
-
- {********** 結果の印字処理 **********}
- procedure Print ;
- var i : range ;
- j : integer ;
- begin
- j := 0;
- write('円周率= ',pai[2]:1,'.');
- for i:=3 to AREA-5 do
- begin
- j := j + 1 ;
- write(pai[i]:1) ;
- if j mod 5 = 0 then write(' ') ; { 5桁ごとに空白 }
- if j=50 then { 50桁ごとに改行 }
- begin
- writeln ;
- write(' ':10) ;
- j:=0
- end
- end ;
- writeln('・・・')
- end ;
-
- {********** 計算処理 **********}
- procedure calc(j1 : integer ; Plus : Boolean) ;
- var j2 : integer ;
- add : Boolean ;
- begin
- p := 0 ;
- j2 := 1 ;
- add := plus ;
- while p < AREA do
- begin
- Divide(term,j1) ; { term := term / j1 }
- Divide(work,j2) ; { work := term / j2 }
- AddSub(add) ; { pai := pai +- work }
- j2 := j2 + 2 ;
- add := not add
- end
- end ;
-
- {********** メイン処理 **********}
- begin {main}
- for i:=0 to AREA do pai[i] := 0 ;
-
- term := pai ; term[1] := 8 ; { 初期値 80.0 }
- calc(25{5*5},true{+-+-・・・・}) ; { マーチンの公式 第1項計算 }
-
- term[0] := 9 ; term[1] := 5 ; term[2] := 6 ;
- for i:=3 to AREA do term[i] := 0 ; { 初期値 956.0 }
- calc(57121{239*239},false{-+-+・・・}) ; { マーチンの公式 第2項計算 }
-
- Print
- end.